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Abstract 

The 1/r Coulomb potential is calculated for a two dimensional system with periodic bound- 
ary conditions. Using polynomial splines in real space and a summation in reciprocal space we 
obtain numerically optimized potentials which allow us efficient calculations of any periodic 
(long-ranged) potential up to high precision. We discuss the parameter space of the optimized 
potential for the periodic Coulomb potential. Compared to the analytic Ewald potential, the 
optimized potentials can reach higher precisions by up to several orders of magnitude. We 
explicitly give simple expressions for fast calculations of the periodic Coulomb potential where 
the summation in reciprocal space is reduced to a few terms. 

1 Introduction 

Most of classical or quantum simulations are using periodic boundary conditions to extrapolate the 
results to the thermodynamic limit of the bulk. Typically, boundary conditions are implemented 
by including replicas of the original system; the potentials are then calculated by considering the 
additional interactions between the particles in the box with all the periodic images of the replicas. 
Whereas for short-range potentials the nearest image convention can often be applied, long-range 
potentials require an additional summation over the Fourier components in reciprocal space due 
to the slow convergence of the contributions of images in real space (see pQ for a recent review of 
how to compute long range potentials within periodic boundary conditions). 

Following Ref. we introduce an optimized periodic potential, represented by two summa- 
tions, one over Fourier and one over real-space components, which can be obtained numerically for 
any potential. We extend the analysis of Ref. [2] to treat two-dimensional systems and concentrate 
particularly on the two-dimensional periodic Coulomb potential. 

For an arbitrary potential v(r), the periodic image potential v pp (r) is defined by summing over 
the interactions between one particle in a box of dimensions L x , L y and the replicas of the other 
particles in periodic space, 



VP 



(r)=5>(r + l) = ]>> k e 4k - r (1) 



Here, 1 are the Bravais lattice vectors {n x L x ,n y L y ) with n x ,n y integers, are the Fourier com- 
ponents of the potential summed over all reciprocal lattice vectors k = 2tt(ti x /L x , n y /L y ) of the 
periodic system. 

For a charge q, the Coulomb potential is given by 

v(r) = - - / dv'y^- , (2) 
r Jv |r-r'| 
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with Fourier components 



^={r ; til w 

in two dimensions where V = L X L V is the volume of the box. A uniform background of opposite 
charge is subtracted to enforce charge-neutrality. Since both summations in real space and in 
reciprocal space of the periodic potential converge slowly, the standard method is to split the 
periodic potentials into two summations, 

Wop (r)=5>(r + l)+ E y^' r W 

1 |k|<K c 

with 

w(r) = for r > R c . (5) 

By definition both summations are converged. In Ref. [2j it has been proposed to use a set of 
basis functions for w(r) and determine numerically their coefficients together with j/k such that 
the difference between the optimized periodic potential v op (r) and the periodic image potential 
v P p(r) is minimized. 

An analytical form for the Coulomb potential was provided long time ago by Ewald using a 
Gaussian charge distribution J2J. It gives 

erfc(ar) 

w (r) = q (6) 



2n crfc(fc/2a) 
aV 



with linv_>o w a (r) — q/r = —q2a/y/n and a is an open parameter which determines the speed of 
convergence in both summations. Due to the exponential convergence of both summations, they 
can be truncated, and R c and K c can be determined to ensure any desired precision. Choosing 
a = y/n/V, both summations roughly converge equally fast, and (R C K C ) is the only parameter 
determining the precision of the truncated Ewald potential. In practice, one typically restricts 
R c < L/2 with L = mm{L x , L y } in order to apply the nearest image convention in real space; the 
precision of the potential then relies on {a, K c }. 

In Ref. |2 it has been shown that a numerical fit of the 3D Coulomb potential reduces consid- 
erably the number of terms in k-space with respect to the analytical Ewald summation in order 
to obtain a comparable precisions. This leads to an important speedup of simulations of charged 
systems. Further, this method is not limited to the Coulomb potential, but can easily be applied 
to any functional form. This is important in ground state quantum Monte Carlo calculations, 
since analytical forms for the potentials have been shown to provide an accurate description of the 
ground state wavefunction of electronic systems however, they typically involve more com- 

plicated (long ranged) functions where no easy analytical break-up can be done. The optimized 
potential of Ref. [2] has the big advantage to be applicable to all type of function; here, we extend 
this method to include two-dimensional systems. 

In the following section, we shortly remind the basic steps necessary to derive the equations of 
the optimized potential which have to be solved numerically, and give the explicit formulas for the 
2D case. The precise numerical evaluation of Bessel functions and their integrals are discussed. In 
section III arc presented the results of the optimized potential. Explicit simple formulas are given 
for the 2D periodic Coulomb potential up to a few percents. 
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2 Method and formulas for 2D 



2.1 General method 

The optimized potential, v opt , is determined by minimizing the absolute error with respect to the 
true periodic potential, wpp, 



x 2 



j2 J a dr l v pp( r ) ~ w °p( r )] 2 • ( 8 ) 
Denoting the Fourier transform of w(r) in the optimized potential, Eq.Q) reads 



v op (r) = J2^ r w k + ^"Sk. (9) 

k |k|<JC e 



and Eq.JSJ is split in two sums: 



x 2 



= E (^k - yk - W k ) 2 + E _ W k) 2 ■ (10) 

|k|<if c |k|>if c 



The first term on the rhs of this equation can be exactly set to zero determining y^, 

Uk = «k - »k for |k| < K c . (11) 

Expanding w(r) — X)<^ c j( r ) using a set of basis functions Cj(r) with Fourier components Cjk, 
Ea. (|ll|l relates the optimal Fourier coefficients ?/k to the optimal coefficients U in real space 



</k 



= Vk~X)*< 5 * for l k l^^- (12) 



The optimal coefficients U can be determined by minimizing the second term of the rhs of Eq. (|10f> 
leading to the following linear equations 

E E C ik C„kiri= E UkCik (13) 
n |k|>AT c |k|>Jf c 

Solving Ea. (|12j) and Ea. (|13(l uniquely determines the optimized potential for any given R c and 
K c . 

2.2 Polynomial Basis set 

Here, we use polynomial splines sitting on a linear grid with m continuous derivatives as basis set 
(assuming a spherical symmetry of the potential). Following Ref.[2], the splines are defined on 
intervals (r^r^+i) with N sp n ne + 1 equally spaced knots starting at the origin and ending at R c , 
Ti — iA with interval A = R c /N sp u ne . The basis functions are Ci a (r) with < a < m are defined 
by imposing 

~O a /30ij. (14) 



drl 3 

The divergence of the potential at the origin is explicitly taken in to acount as follow 



Cia{r) 

i=0 a=0 



(r)= E E'^ 



where C = 1 for the Coulomb potential and C — for regular potential at the origin. Thus, the 
basis functions Ci a {r) are piecewise polynomial of order 2m + 1 



c (r) ,A«E:TSo„(¥)" : n<r<r l+1 
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and zero for |r — 7"j| > A. The constraints at r — fj fixe half of the S-elements 

S an — — r<Wi) for 0<a,rt<m. (17) 
m 

The constraints at r = Vi±\ gives 

" 1 

S a , n+m +i =-J2 i M ~ l ) k n ( _ fc v for 0<a,n<m (18) 
where M _1 is the inverse of the quadratic matrix 

M ak = / ■ 7 / j M| for 0<a,fe<m. (19) 
(m + 1 + a — k)\ 

The required Fourier coefficients 5i a k are given by 

2m+l 

Cfa* = A" ^ S an (D+ n + (-ir +n D7 kn ) (20) 



where 



1 



71=0 



^ = ^1 dre " ikrr " C (, (21) 
= ±17XTrEG")(-^7 dr^e-** (22) 



VA" 

9_ n rri±i 

= ± vM^®(-n) n - j drr^- c J (kr) (23) 

with Jo(x) is the Bessel function of zero order and (™) denotes the binomial coefficients. The 
moments of Jo can be obtained from its first two moments and a reccurence relation: 

dxJ (x) = xJ (x) + ^(H {x)J 1 (x)-H 1 (x)J (x)) (24) 
dxx Jq(x) = xJi(x) (25) 
dxx n J (x) = x n Ji{x) + (n- l^^Mx) -{n- I) 2 [ dxx n - 2 J (x). (26) 



Here, Ji is the Bessel function of first order and Hi (H2) is Struve's function of first (second) 
order. However, contrary to the 3D case, we have not found any "machine precision" routine to 
compute the Struve's function or the integral of Jo due to the oscillatory behavior of the integrand. 
In Appendix A, we briefly describe how to evaluate "precisely" the integral of Jq. 



3 Results for the 2D periodic Coulomb potential 

In this section, we present the results of the optimized potential for the Coulomb 1/r potential 
in two dimensions, Eq.J2J) and Eq.Q in a square box of length L = L x = L y . Hcrmite splines 
of fifth order represent the real space part of the optimized potential insuring two continuous 
derivatives (m — 2 in Ea llfijl . The 1/r divergence at the origin is accounted for by setting C = 1 
with the constraint ti = o,a=o = Q] symmetry further imposes the absence of any linear term linear: 
U=o. a =2 = 0. Both constraints can be easily included, by solving the linear equations 

53 A-iajptjp =bia, < a < m nderv , < i < N sp u ne (27) 

0\/3)#{(0,0),(0,2)} 
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for (i,a) ^ {(0,0), (0,2)} with 



K m K m 
Aiu,j[i = ^ CiakCjpk ha — ^ (Vfc — CoofetoO ) Ci a k (28) 
k=K c k=K c 

instead of Ea. i|13|) . Appendix B describes how to impose the Madelung constant of the lattice 
through the constant term toi ■ 

The number of spline intervals, N sp u ne , and the cutt-off K c in the reciprocal space are the 
open parameters of the model. In Ea. l|28|l . K m has to be large enough to insure that each sum 
is converged. Since it is an important parameter which determines the conditioning of the linear 
equation we first discuss its influence on the stability of the solution. 

The results of the optimized potential are compared with the "exact" periodic Coulomb po- 
tential obtained from the Ewald formula including a summation over many images in real space. 
Thus machine precision is easily reached for this reference potential. Values of the potential are 
given in units of q/L, they are independent of the size of the box. 



3.1 Extrapolation of K m 

Results are very sensitive to K m which determines the convergence of the matrix elements of A 
and b, Ea. l|28() . From Ea. (|20l) and Eci. (|23[) . one finds the dominant behavior in the limit of k — + oo: 
toocook — Vk ~ 0(k~ 7 / 2 ) and Ci a k ~ 0(/c~ 7 / 2 ) for (i, a) ^ (0,0). Thus the matrix elements in A 
and 6, Eq.J^HJ, are of order fc~ 7 . In the large k limit, the truncation therefore introduces an error 
of order K~ 5 . 

Considering the asymptotic expansions of ci a k for large k, one might be able to extend the 
summation analytically to infinity using the leading order terms. However, in 2D, the difference 
between the discrete summation over reciprocal lattice vectors and the continuous integration is 
comparable to the correction of the continuous integral. Therefore, we find no improvement by 
adding those corrections. Thus, contrary to the 3D case, no analytical continuations are used here. 

Next, we consider the stability of the solution with respect to variations of K m . Since the short 
wavelength cut-off destroys the information on small distance behavior of the potential, we expect 
that there is a maximum number of splines N™f^ e , after which the solution of the linear equation 
will become unstable. Roughly, the resolution in real space is limited by AK m / (2m + 2) >^ 2ir, 
and we get a maximum number of spline knots N™^ ne + 1 with 

N max < K m R c 
^ spline S -7—, — TT 

F 47r(m + 1) 

Using less terms in the summation in reciprocal space, the matrix Ai a jp in Ea. (|28|l becomes 
ill conditioned. The conditioning of the linear system can be estimated by comparing the norm of 
the obtained solution t t , \\t\\oo = max{|tj|,i} with HsH^ = max{| Y,j Aijtj-h], i}. If | |s| |oo/| |t| \ x 
is of order one, the system is dominated by numerical round-off errors. This indicates that cither 
the value of K rn is too small or the number of splines is too large so that no improvement can be 
reached by increasing N sp u ne further. 



3.2 Accuracy of the optimized potential 

We now study the accuracy of the optimized potential in the N sp u ne - K 2 -plane for fixed R c /L = 
0.5 (nearest image convention). The \px i ~ in units of q/L corresponds to the average error and 
is shown in Fig. ^ We see that there is an optimum line in the range of parameters considered. 
We also note that there is a difficulty to decrease the precision below 10~ 10 , mainly because the 
linear system, Eq.JSSJ) becomes more and more ill-conditioned as N sp u ne and K 2 increase. 

Figure |21 shows the difference to the true periodic potential for N sp n ne = 30 varying K 2 
and compares it with the values of the best nearest image potential using the analytic Ewald 
expressions, Eq.Q with a = K c /L. For the range of interest in order to speed up simulations, 
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Figure 1: Contour plot of the mean error of the optimized potential versus the number of shells 
N% (N c = K c L/2ir) and the number of splines N sp u nes (see definition of Eq.Q). The maximum 
distance in r-space is chosen to satisfy the nearest image convention in simulations, R c /L = 0.5. 



the optimized potential is always better by at least one order of magnitude. However, if very high 
precision is needed, better than 10 -10 , the Ewald method is preferable to the present optimized 
procedure, even if it is much more cpu-time consuming, since it has no stability problems. 

The number of splines increases the precision of the potential inside the circle of radius R c . 
Around the corner of the box, only the number of fc-shells can improve the optimized potential. At 
fixed number of shells, one reaches rapidly an optimum number of splines after which increasing 
the number of splines has no more effect. This is seen in Fig|2]by the straight vertical lines. At 
fixed number of splines, increasing the number of fc-shells first improves strongly the optimized 
potential, but later the improvement almost saturates. Thus the optimum choice is to take the 
parameters roughly along the diagonal in FigEI 

Note that an intermediate precision of 10~ 3 to 10~ 6 is reached with very small values of the 
optimized potential parameters. Therefore simple analytical expressions allow us fast evaluations 
of the periodic Coulomb potential with intermediate precisions. 



3.3 Simple expressions for intermediate precision 

If we include only the first five wavevectors with |k| < K c = 2ir/L and use N sp u ne = 2, the 
minimum to obtain a smooth curve going to zero at half of the box size, a mean precision of 1% 
is obtained, and a maximum error of around 2% (see Table l3~3*|) . Increasing the number of splines 
improves slightly the precision. However, extending for N S pUne = 2 the reciprocal summation 
using K c — in/L, the precision decreases to 0.1%, with a maximum error of 0.4% (see Table l3~3*|) . 
The short-range part of the optimized Coulomb potential writes 

w[r) -l{ Eto . : 0<r<L/4 

whereas the long-range part is given by 

y( r ) = ^ {Vo + 2yi [cos(x) + cos(y)] + 4y 2 cos(x) cos(y) + 2y 4 [cos(2x) + cos(2y)]} (31) 

where x — 2irx/L, y — 2iry/L. See Fig|3]for a comparison of these simple expression with the 
"exact" periodic potential. 
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Figure 2: Comparison of the precision between the standard Ewald method (with a — K c /L) 
(triangle) and the the optimized potential versus N c = K c L/2tt (square). The average error \ 1S 
given in units of q/L. The maximum distance in r-space is chosen to satisfy the nearest image 
convention, R c /L = 0.5. 
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Figure 3: Difference between the simple expression of the optimized potential, v op , with the exact 
Coulomb potential, v pp ] dooted lined stands for v op using K c = 47r/L, full line for K c = 8tt/L. O 
stands for the origin, A for the middle of the square side and B for the corner of the square. 
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5 




1 


-0.819506 





0.169304 


-0.146967 


0.0777952 




0.280626 


-0.510485 


0.404063 


-0.955541 


1.3377 


-0.556365 


n 2 





1 










j/n 2 


-1.11863 


0.124098 











Table 1: Optimized Potential parameters, Ea. H30l and Ea. l31fl . Top : short range real space 
parameters for N sp u ne = 2. Bottom : reciprocal space parameters where n = kL/2it. The mean 
precision is about 2%. 
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-1.09583 





0.30778 


-0.0359887 


-0.0266302 


a> 


0.149336 


-0.449592 


0.441105 


-0.119121 


-0.0333851 


0.0116568 


n 2 





1 


2 


4 






fin* 


-0.870938 


0.262177 


0.0715766 


0.00474028 







Table 2: Same as Table EOl The mean precision of this model is about 0.1%. 



Both explicit expressions are roughly one order of magnitude better than the corresponding 
Ewald potentials. Even if not extremely precise, the expressions should extrapolate much better 
to the thermodynamic limit than any truncated potential using only nearest image convention 
in real space. Further the real space part of the optimized potential vanishes at R c = L/2 by 
construction without introducing any discontinuity in the potential and the derivatives at this 
point. 

4 Conclusion 

We have shown that for the two dimensional Coulomb potential, the numerically optimized po- 
tential can obtain a much higher precision compared to the analytical Ewald potential summing 
over the same number of terms in reciprocal space. Therefore, the computational effort for many- 
body simulations involving long range potentials can be significantly reduced using an optimized 
potential. 

For a pair potential the computational cost to evaluate the real space contribution to the total 
potential is ~ NN c /2 where N is the total number of particle, iV c = irR 2 p is the "number of close 
neighbors" and p = N/V the mean particle density. Since the number of k— vectors increases as 
the volume in reciprocal space, the cost of the Fourier summation is roughly ~ NttK 2 V/4:Tt 2 , and 
there is an optimum value for each system size, which scales as R c ~ K~ x ~ TV 1 / 4 ~ L 1 / 2 in the 
limit of a large particle number, so that the computational cost in reciprocal space ~ N 3 / 2 roughly 
equals that in real space. We further note, that in the limit of a large system, the total cost for 
evaluating the Coulomb potential using real and reciprocal space summations is always favorable 
compared to any truncated potential with minimum image convention, which scales as ~ iV 2 . 

Apart from a potential speed-up of simulations involving charged particles, the big advantage 
of the optimized potential is its flexibility to split-up any (long-ranged) function into a real-space 
and a reciprocal space contribution. 

Appendix A 

In this Appendix, we describe how to evaluate precisely the integral of the Bessel function Jo (ir) . 
Beginning with the evaluation of the Bessel function J$, three domains are defined. Around the 
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origin, < x < x\ the Bessel function is accurately evaluated form the absolutely convergent 
series representation [H] 

With standard double precision, a precision of 10~ 14 is obtained up to x\ — 10 by summing all 
terms whose absolute value is larger than 10~ 18 . For large arguments, X2 < x < 00, the asymptotic 
expansion is used: jH] 

J (x) = \ —{P(x)cos(x - ir/4) + Q(x)sin(x - n/4)} , x 2 < x < 00 (33) 

V TTX 

with 

P (T) = 1 + y-(^ 2 fr 1 ( 2 - + 1 ) 2 Q, X) = T tH fr ( 2 - + D 2 (34) 

As usual, the summation is stopped when the absolute value of the running term of the series 
starts to increase. With the asked precision of 10 -14 , one finds x-i — 16. In the intermediate 
region x\ < x < X2, the Chebyshev-pade approximant are calculated using Maple. This strategy 
can be extended to any desired precision by cutting this interval in pieces. 

The series expressions for Jq(x) are then analytically integrated to calculate J dxJo(x). Un- 
fortunately, the asymptotic expression is less converging, giving X2 = 30 for a precision of 10~ 14 . 
The Chebyshev-pade approximant in the interval [10, 30] is calculated at order 36, thanks to Maple: 
with(orthopoly): with(numapprox): Digits:=20; IJ:=int(BesselJ(0,x),x); IJCh:=eval(chebpade(IJ, 
x=10.. 30,36)): convert(suhs(x=10*(X+2) y IJCh),horner); where X = (x-20)/10. The gsl routine 
have been used for Jq and J\ [7] . 



Appendix B 

It is also possible to fix the constant term at the origin in the optimized Coulomb potential, to\ in 
Ea. (|15fl . by imposing the Madelung constant of the underlying lattice, VMad = Hm r _>o ( v pp( r ) — <z/ r ) 
Using the Ewald expressions, Eq.J7|), the Madelung constant writes 

v Mad = £ ™ Q (l)+E^- 2 \/! ( 35 ) 

1/(0,0) k V 

which can be calculated with high precision choosing a — i/tt/L; for the square lattice VMad = 
— 3.90026492000195g/i. For the optimized potential, imposing the 1/r divergency with C = 1, 
the Madelung verifies VMad = £01 + J2k<K V k - Since yu is coupled to all U a by Ea. (|12Jl . this 
constraint leads to modifications of Ea. l|28fl . 

^ Aiajptjp = b ia , < a < m, < i < N spUne (36) 

(j,«/{(0,0),(0,1),(0,2)} 

for (i, a) ^{(0,0), (0,1), (0,2)} with 

K m 

A-ia.jfl = ^ GiakGifik (37) 



k=K a 
K m 



bia — ^2 ^fc ~ CQOktoO — ClOfc -. _ Mad ^2 tjp(Cjl3q — Cjf3q) \ Ciak (38) 

k=K c \ 1 ~~ ^q<K c C Wq jm0tl) J 

X k = X fe +c 10fc l ^ K " X ! (39) 

1 _ l^q<K c C Wq 
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instead of Ea. (|13[) . Here X denotes either Ci a or v. Including the Madelung term as a constraint 
improves slightly the solution. 
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